;;########################################################################
;; nway.lsp
;; Copyright (c) 1993-99 by Forrest W. Young
;; Nway Analysis of Variance object. Can fit main effects or two-way 
;; interaction models to n-way balanced or unbalanced univariate data.
;; Data must be complete (no empty cells)
;;
;; Uses Tierney's regression proto with Scheffe dummy variables.
;;########################################################################

;;constructor function

;fwy 4.28+ changed to allow analysis of unbalanced data
(defun nway-model 
  (data &key (print t) way-names class-names interaction)
"Args: (data &key (print t) :group-names :way-names :interaction)
DATA: List of lists of compound-data, there being one sub-list for every cell of the design.  No empty, missing or nil cells permitted.  Unbalanced data (unequal cell frequencies) are permitted. 
WAY-NAMES:   List of n character strings naming the ways of the design.  
CLASS-NAMES: List of n lists of character strings, one for each of the n ways of the design.  Each list contains the names of the classes (levels) of a way.
INTERACTION: Nil for no interaction terms, t for two-way interaction terms. Higher-way interactions not computed. CLASS-NAMES and WAY-NAMES must be specified. INTERACTION is optional and cannot be used with one-way data or un-replicated two-way data."
  (let ((data (mapcar #'(lambda (x) (coerce x 'list)) data))
        (m (send nway-model-proto :new
                 data way-names class-names interaction)))
    (if print (send m :display))
    m))

(defproto nway-model-proto 
  '(data nobs nways nclasses    cellfreqs interaction
         way-names  class-names source-names
         source-sum-of-squares  source-degrees-of-freedom
         design-matrices        indicator-matrices 
         main-effects-matrices  interaction-matrices) 
  '() regression-model-proto)

(defmeth nway-model-proto :isnew 
  (data way-names class-names interaction)
  (send self :data data)
  (send self :nobs (length (combine (send self :data))))
  (send self :cellfreqs (mapcar #'length (send self :data)))
  (send self :class-names class-names)
  (send self :nclasses (mapcar #'length (send self :class-names)))
  (send self :way-names way-names)
  (send self :nways (length (send self :way-names)))
  (send self :make-main-effects-matrices)
  (send self :main-effects-matrices (send self :design-matrices))
  (when (= (send self :nways) 1) (setf interaction nil))
  (when interaction 
        (send self :interaction t)
        
        (send self :make-interaction-matrices)
        (send self :design-matrices 
              (concatenate 'list (send self :design-matrices)
                           (send self :interaction-matrices))))
  (send self :y (combine data))
  (send self :x (apply #'bind-columns (send self :design-matrices)))
  (send self :make-source-names)
  (send self :make-predictor-names)
  (call-next-method)
  (send self :intercept t)
  (send self :make-ssq-df)
  )

(defmeth nway-model-proto :nways (&optional (number nil set))
"Args: (&optional number)
Sets or returns number of ways."
  (if set (setf (slot-value 'nways) number))
  (slot-value 'nways))

(defmeth nway-model-proto :nobs (&optional (number nil set))
"Args: (&optional number)
Sets or returns number of observations."
  (if set (setf (slot-value 'nobs) number))
  (slot-value 'nobs))

(defmeth nway-model-proto :nclasses (&optional (number-list nil set))
"Args: (&optional number-list)
Sets or returns a list of the number of classes for each way."
  (if set (setf (slot-value 'nclasses) number-list))
  (slot-value 'nclasses))

(defmeth nway-model-proto :cellfreqs (&optional (number-list nil set))
"Args: (&optional number-list)
Sets or returns a list of the frequency of each cell."
  (if set (setf (slot-value 'cellfreqs) number-list))
  (slot-value 'cellfreqs))

(defmeth nway-model-proto :way-names (&optional (string-list nil set))
"Args: (&optional string-list)
Sets or returns way-names."
  (if set (setf (slot-value 'way-names) string-list))
  (slot-value 'way-names))

(defmeth nway-model-proto :class-names (&optional (names nil set))
"Method args: (&optional names)
Sets or returns class names."
  (if set (setf (slot-value 'class-names) names))
  (slot-value 'class-names))

(defmeth nway-model-proto :interaction (&optional (logical nil set))
"Method args: (&optional logical)
Sets or returns t or nil indicating whether two-way interaction terms used."
  (if set (setf (slot-value 'interaction) logical))
  (slot-value 'interaction))

(defmeth nway-model-proto :design-matrices (&optional (matrix-list nil set))
"Method args: (&optional matrix-list)
Sets or returns the design matrices."
  (if set (setf (slot-value 'design-matrices) matrix-list))
  (slot-value 'design-matrices))

(defmeth nway-model-proto :indicator-matrices 
  (&optional (matrix-list nil set))
"Method args: (&optional matrix-list)
Sets or returns the indicator matrices."
  (if set (setf (slot-value 'indicator-matrices) matrix-list))
  (slot-value 'indicator-matrices))

(defmeth nway-model-proto :main-effects-matrices 
  (&optional (matrix-list nil set))
"Method args: (&optional matrix-list)
Sets or returns the main effects matrices."
  (if set (setf (slot-value 'main-effects-matrices) matrix-list))
  (slot-value 'main-effects-matrices))

(defmeth nway-model-proto :interaction-matrices 
  (&optional (matrix-list nil set))
"Method args: (&optional matrix-list)
Sets or returns the interaction matrices."
  (if set (setf (slot-value 'interaction-matrices) matrix-list))
  (slot-value 'interaction-matrices))

(defmeth nway-model-proto :data (&optional (list nil set))
"Message args: (&optional list)
Sets or returns the data."
  (if set (setf (slot-value 'data) list))
  (slot-value 'data))

(defmeth nway-model-proto :source-names (&optional (string-list nil set))
"Message args: (&optional string-list)
Sets or returns the source names."
  (if set (setf (slot-value 'source-names) string-list))
  (slot-value 'source-names))

(defmeth nway-model-proto :source-sum-of-squares 
  (&optional (number-list nil set))
"Args: (&optional sum-of-squares)
Sets or returns a list of the sum-of-squares for each source."
  (if set (setf (slot-value 'source-sum-of-squares) number-list))
  (slot-value 'source-sum-of-squares))

(defmeth nway-model-proto :source-degrees-of-freedom 
  (&optional (number-list nil set))
"Args: (&optional number-list)
Sets or returns a list of the degrees-of-freedom for each source."
  (if set (setf (slot-value 'source-degrees-of-freedom) number-list))
  (slot-value 'source-degrees-of-freedom))

(defmeth nway-model-proto :make-source-names ()
"Method args: none
Makes source names from way names according to ANOVA specifications."
  (let* ((source-names (send self :way-names))
         (nways (send self :nways))
         (knt nways)
         )
    (when (send self :interaction)
          (setf source-names 
                (combine source-names 
                         (repeat " " (/ (* nways (- nways 1)) 2))))
          (dotimes (i (- nways 1))
                   (dolist (j (iseq (+ i 1) (- nways 1)))
                           (setf (select source-names knt)
                                 (strcat (select source-names i) "*"
                                         (select source-names j)))
                           (setf knt (+ knt 1)))))
    (dotimes (i (length source-names))
             (if (> (length (select source-names i)) 20)
                 (setf (select source-names i) 
                       (subseq (select source-names i) 0 20))))
    (send self :source-names source-names)))
 
(defmeth nway-model-proto :make-predictor-names ()
"Method args: none
Makes predictor names from class names according to ANOVA specifications."
  (let ((class-names (send self :class-names))
        (way-names (send self :way-names))
        (nways (send self :nways))
        (nclasses (send self :nclasses))
        (i-classes nil)
        (last-i-class nil)
        (i-class nil)
        (j-class nil)
        (pred-names (send self :predictor-names))
        (k 0)
        (l 0)
        (m 0)
        (kinc 0)
        (linc 0))
    (dotimes (i nways)
             (setf i-classes (select class-names i))
             (setf last-i-class (subseq (car (last i-classes)) 0 1))
             (dotimes (j (- (select nclasses i) 1))
                      (setf j-class (subseq (select i-classes j) 0 1))
                      (setf (select pred-names k)
                            (strcat (select way-names i)
                                    "[" j-class "-" last-i-class "]"))
                      (setf k (+ 1 k))))
    (setf m k)
    (when (send self :interaction)
          (dotimes (i (- nways 1))
             (setf kinc 0)
             (dotimes (kp i)
                (setf kinc (+ kinc (select nclasses kp) -1)))
             (setf i-class (select nclasses i))
             (dolist (j (iseq (+ i 1) (- nways 1)))
                (setf j-class (select nclasses j))
                (setf linc 0)
                (dotimes (lp j)
                   (setf linc (+ linc (select nclasses lp) -1)))
             (dotimes (k (- i-class 1))
                (dotimes (L (- j-class 1))
                   (setf (select pred-names m) 
                         (strcat (select pred-names (+ k kinc)) "*"
                                 (select pred-names (+ l linc))))
          (setf m (+ m 1)))))))
    (send self :predictor-names pred-names)))

(defmeth nway-model-proto :make-main-effects-matrices ()
"Computes indicator and design matrices for the main effects of balanced or unbalanced tabular data having any number of ways. Each matrix has one row for each observation. The indicator matrices have a column for each level but are not full rank. Their elements are 1's and 0's, with a 1 indicating which level the observation belongs in.  The design matrices have one column less than the number of levels and are full rank. There elements are 1's, 0's and -1's, as suggested by Scheffe."
  (let* ((nways (send self :nways))
         (nobs  (send self :nobs))
         (freqs nobs)
         (nlevels (send self :nclasses))
         (cellfreqs (send self :cellfreqs))
         (nthislevel nil)
         (nproduct 1)
         (nrepeats 1)
         (thisobs 0)
         (thiscell 0)
         (thisfreq 0)
         (design nil))
    (dotimes (i nways)
             (setf nproduct (* nproduct (select nlevels i))))
    (dotimes (i nways)
;determine cell frequencies for way i of the design
             (setf nthislevel (select nlevels i))
             (setf nproduct (/ nproduct nthislevel))
             (setf freqs (repeat 0 (* nthislevel nrepeats)))
             (dotimes (ii (* nthislevel nrepeats))
                      (dotimes (jj nproduct)
                               (setf (select freqs ii) 
                                     (+ (select freqs ii)
                                        (select cellfreqs thiscell)))
                               (setf thiscell (+ thiscell 1))))
             (setf thiscell 0)
;create reduced-rank indicator matrix for way i of the design
             (setf design 
                   (make-array (list nobs nthislevel) :initial-element 0))
             (dotimes (k nrepeats)
                      (dotimes (j nthislevel)
                               (dotimes (l (select freqs thisfreq))
                                        (setf (aref design thisobs j) 1)
                                        (setf thisobs (+ thisobs 1)))
                               (setf thisfreq (+ thisfreq 1))))
             (setf nrepeats (* nrepeats nthislevel))
             (send self :indicator-matrices 
                   (make-matrix-list (send self :indicator-matrices) design))
             (setf thisobs 0)
             (setf thisfreq 0)
;create full-rank design matrix for way i of the design
             (setf design (select (- design (transpose 
                          (matrix (list nthislevel nobs ) 
                                  (repeat (col design (- nthislevel 1)) 
                                          nthislevel)))) 
                                  (iseq nobs) (iseq (- nthislevel 1))))
             (send self :design-matrices 
                   (make-matrix-list (send self :design-matrices) design)))))

(defmeth nway-model-proto :make-interaction-matrices ()
"Args: none
Computes and stores in object slot all two-way interaction design matrices."
  (let ((n (length (send self :way-names))))
    (when (> n 1)
          (dolist (i (iseq (- n 1)))
                   (dolist (j (iseq (+ i 1) (- n 1)))
                           (send self :interaction-matrices
                              (make-matrix-list 
                                (send self :interaction-matrices)
                                (send self :twoway-interaction-matrix j i)))
                           )))))

(defmeth nway-model-proto :twoway-interaction-matrix (wayi wayj)
"Args: WAYI WAYJ
WAYI and WAYJ are integers specifying a way of the design.  Computes and returns the two-way interaction design matrix for the two ways."
  (let ((designi  (select (send self :design-matrices) wayi))
        (designj  (select (send self :design-matrices) wayj))
        (classesi (select (send self :nclasses) wayi))
        (classesj (select (send self :nclasses) wayj))
        (interact nil)
        (nobs     (send self :nobs))
        (interaction nil))
    (dotimes (i (- classesj 1))
             (setf interaction (make-matrix-list interaction (transpose 
                   (matrix (list (- classesi 1) nobs)
                           (map-elements #'* 
                                   (repeat (col designj i) (- classesi 1)) 
                                   (transpose designi)))))))
    (apply #'bind-columns interaction)))

(defmeth nway-model-proto :make-ssq-df ()
"Message args: none
Computes and saves the sum of squares and degrees of freedom of each source in the anova for balanced or unbalanced data. SSQ and DF are saved as lists with the overall model values first, the sources next, with the error term last."
  (let* ((source-df (combine (mapcar #'size (send self :design-matrices))))
         (source-df 
          (select source-df (+ 1 (* 2 (iseq (/ (length source-df) 2))))))
         (cumdf (cumsum source-df))
         (b (select (send self :coef-estimates) (iseq 1 (first source-df))))
         (x (first (send self :design-matrices)))
        ;(means (matmult x b))
        ;(sumsq (sum (^ means 2)))
         (indices nil)
         (m1 (select (send self :xtxinv) 
                     (iseq 1 (first source-df)) (iseq 1 (first source-df))))
         (sumsq-new (matmult b (inverse m1) b)
          (/ (^ b 2) (aref (send self :xtxinv) 1 1)));fwy 4.28+
         (modeldf (- (length (send self :coef-estimates)) 1))
         (errordf (send self :df))
         (errorss (send self :sum-of-squares))
         (modelss (- (send self :total-sum-of-squares) errorss))
         )
    (send self :source-sum-of-squares (list modelss sumsq-new))
    (send self :source-degrees-of-freedom 
          (combine modeldf source-df errordf))
    (when (> (length source-df) 1)
          (dolist (i (iseq 1 (- (length source-df) 1)))
                  (setf b (select (send self :coef-estimates) 
                                  (iseq (+ (select cumdf (- i 1)) 1)
                                        (+ (select cumdf (- i 1))
                                           (select source-df i)))))
                  (setf x (select (send self :design-matrices) i))
                ; (setf means (matmult x b))
                ; (setf sumsq (sum (^ means 2)))
                  (setf indices (iseq (1+ (select (cumsum source-df) (1- i)))
                                      (select (cumsum source-df) i)))
                  (setf m1 (select (send self :xtxinv) indices indices))
                  (setf sumsq-new (matmult b (inverse m1) b))
                  (send self :source-sum-of-squares 
                        (combine (send self :source-sum-of-squares) sumsq-new))))
    (send self :source-sum-of-squares 
          (combine (send self :source-sum-of-squares) errorss))))

(defmeth nway-model-proto :ls-means ()
  (let* ((grouped-data-and-labels 
         (send model-obj :grouped-data (1- sel) point-labels))
        (grouped-data (first grouped-data-and-labels))
        (permuted-point-labels 
         (combine (second grouped-data-and-labels)))
        (means (mapcar #'mean grouped-data))
        )
    ))

(defmeth nway-model-proto :display 
  (&optional w (details t) model equation interpret (flush t))
"Message args: (w details model)
Prints the standard regression summary plus the effects table in window W.  If W is nil, prints it on standard output (listener)."
  (let* ((source-df (send self :source-degrees-of-freedom))
         (source-ss (send self :source-sum-of-squares))
         (error-df  (first (last source-df)))
         (error-ss  (first (last source-ss)))
         (total-df  (sum (rest source-df)))
         (total-ss  (sum (rest source-ss)))
         (model-df (- (length (send self :coef-estimates)) 1))
         (model-ss (- total-ss error-ss))
         (r-squared (send self :r-squared))
         (nways (send model :nways))
         (source-result)
         (results)
         (titleline "~%Source             Sum-of-Squares   df  Mean-Square    F-Ratio   P-Value    R-Square~%")
         )
    (if (and (not interpret) details)
        (call-next-method w (member 1 details) (member 2 details) 
                          (and (member 3 details) (> nways 1)))
        (display-string (format nil "~%") w))
    (unless interpret
            (display-string (format nil "~2%ANALYSIS OF VARIANCE: EFFECTS TESTS                            Signficance  Strength") w)
            (display-string (format nil titleline) w))
    (setf source-result (send self :anova-table-line w 
          (first (send self :source-names))
          (second source-ss) (second source-df) 
          error-ss error-df total-ss interpret))
    (setf source-result (combine source-result (first (send self :source-names))))
    (when (> (length source-df) 3)  
          (dolist (i (iseq 1 (- (length source-df) 3)))
                  (send self :anova-table-line w 
                        (select (send self :source-names) i)
                        (select source-ss (1+ i))
                        (select source-df (1+ i)) 
                        error-ss error-df total-ss interpret))
          (setf results (send self :anova-table-line w 
                              "All Sources" model-ss model-df 
                              error-ss error-df total-ss
                              interpret)))
    (when (not results) (setf results source-result))    
    (unless interpret
            (send self :anova-table-line w "Error" error-ss error-df)
            (send self :anova-table-line w "Total" total-ss total-df)
            (setf unique-ss (- (send self :total-sum-of-squares) total-ss))
            (when (and (> (length (send self :source-names)) 1) (> unique-ss 0))
                  (send self :anova-table-line w "Unique" unique-ss))
           ; (when (not details) 
           ;       (display-string (format nil "~%OVERALL STRENGTH: ") W)
           ;       (display-string (format nil "R-Squared             = ~4,2f~%" 
           ;                               r-squared ) w))
            )
    (when interpret
          (send self :interpret-nway-model model results source-result 
                model-df error-df (second source-df) nways flush))))

(defmeth nway-model-proto :interpret-nway-model 
  (model results source-result model-df error-df source-df nways &optional (flush t))
  (let* ((p-val-str (strcat "p="(second results)))
         (p-src-val-str (strcat "p="(second source-result)))
         )
    (when flush (interpretation-window (send model :title) :flush flush))
    (if (> nways 1)
    (paste-interpretation (format nil "F-TEST - Full Model:~%For these data we performed an F-Test having ~d and ~d degrees of freedom. This test resulted in F=~3,2f, ~a. This means that when the null hypothesis of no effect for any source in the model is true, the probability is ~a of obtaining a sample of data with an effect that is as extreme or is more extreme. Thus, of many repetitions of this sampling design, the probability is ~a that we would be wrong when we reject the null hypothesis of no effect in favor of the alternative hypothesis that at least one pair of treatment means is different."
        model-df error-df (first results) p-val-str p-val-str p-val-str) flush))

    (paste-interpretation (format nil "~2%F-TEST") )
    (if (> nways 1) (paste-interpretation (format nil " - Each Model Source:")))
    (paste-interpretation (format nil "~%For ~a we performed an F-Test having ~d and ~d degrees of freedom. This test resulted in F=~3,2f, ~a. This means that when the null hypothesis of no effect for ~a is true, then the probability is ~a of obtaining a sample of data with an effect that is as extreme or is more extreme. Thus, over many repetitions of this sampling design, the probability is ~a that we would be wrong when we reject the null hypothesis of no effect in favor of the alternative hypothesis that at least one pair of ~s means is different."
              (fourth source-result) source-df error-df 
              (first source-result)  p-src-val-str
              (fourth source-result) p-src-val-str p-src-val-str 
                                    (fourth source-result)) )
    (when (> nways 1) 
          (paste-interpretation (format nil " A similar interpretation can be made for other model sources.") )
          
          (paste-interpretation (format nil "~2%STRENGTH OF RELATIONSHIP - Full Model:~%Based on this sample of data, we estimate that ~2d% of the variability in ~a is related to variability in all sources. This estimate would vary over repeated samples of data."
               (fuzz (* 100 (third results)) 0) (send model :response)) ))
  
    (paste-interpretation (format nil "~2%STRENGTH OF RELATIONSHIP") )
    (if (> nways 1) (paste-interpretation (format nil " - Each Model Source") ))
    (paste-interpretation (format nil "~%Based on this sample of data, we estimate that ~2d% of the variability in ~a is related to variability in ~a. This estimate would vary over repeated samples of data."
                  (fuzz (* 100 (third source-result)) 0) (send model :response) 
                                  (fourth source-result)) )
    (if (> nways 1) 
        (paste-interpretation (format nil " A similar interpretation can be made for other model sources.")) )
    (when flush (show-interpretation))
    ))